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I O , Abstract. - The Poisson-Boltzmann approach gives asymptotically exact counter-ion density 

profiles around charged objects in the weak-coupling limit of low valency and high temperature. 
In this paper we derive, using field-theoretic methods, a theory which becomes exact in the 
K^ ^ opposite limit of strong coupling. Formally, it corresponds to a standard virial expansion. 

I i»», , Long-range divergences, which render the virial expansion intractable for homogeneous bulk 

l^^_ ■ systems, are shown to be renormalizable for the case of inhomogeneous distribution functions 

_«>. ' by a systematic expansion in inverse powers of the coupling parameter. For a planar charged 

f^y. I wall, our analytical results compare quantitatively with extensive Monte-Carlo simulations. 
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Recent years have witnessed a revival of the interest in classical charged system s [|1[, |2[] . Spe- 
cific attention has been paid to the failure of the Poisson-Boltzmann (PB) approach^, y, y, g, 0, 
which is known to give reliable results only in the limit of low- valency ions or high temperatures. 
pH ' Corrections to PB have been attributed to correlations between ions, or, more precisely, corre- 

O , lated ion-density fluctuations, and, if present, additional non-electrostatic interactions between 

O ' ions. These corrections are particularly important for the interaction between macroscopic 

similarly charged objects, where they can lead to attractions |^, ^, ||, H, ||. In as much as 
the PB approach is accurate for weakly charged systems, no systematic theory is available for 
the distribution of counter-ions around charged objects in the opposite limit of high-valency 



C^ ' ions; moreover, it was not clear whether such a limit exists and whether it is physically 

meaningful. In this paper we show, using field-theoretic methods, that while PB corresponds to 
the asymptotically exact theory in the weak- coupling limit, (corresponding to low- valency ions 
or high temperatures), our novel strong- coupling theory becomes asymptotically exact in the 
opposite limit of high- valency ions or low temperatures and constitutes a physically sound limit. 
For the case of a planar charged wall, we give explicit results for the asymptotic density profile 
in the strong-coupling limit. We also have performed extensive Monte-Carlo (MC) simulations 
of this system. The resulting density profiles agree for weak and for strong coupling with 
predictions from PB theory and our strong-coupling theory, respectively. The strong-coupling 
limit is experimentally easily reached at room temperatures with highly charged walls and/or 
multivalent counter ions and thus relevant from the application point of view. 
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To proceed, consider the Hamiltonian of a system of N ions of valency q at an impenetrable 
and oppositely charged wall of number density of surface charges CTj, 

j<k' ^ "' 3=1 

where £b = e^ / AirekBT is the Bjerrum length which measures the distance at which two unit 
charges interact with thermal energy; in water £3 ~ 0.7nm. The dielectric constant is assumed 
to be homogeneous throughout the system. The Gouy-Chapman length, fj, = l/2'!Tqas£B, 
measures the distance from the wall at which the potential energy of an isolated ion reaches 
thermal energy. Rescaling all lengths by /x according to r = iif, the Hamiltonian reads 

j<k ■' 3 = 1 

and thus only depends on the coupling parameter S — 2Trq^i^as. Using the fact that the 
typical distance of an ion from the wall is in reduced units z ^ 1 (which holds both in the weak 
and strong-coupling limits, as we will show below) the confinement energy, the second term in 
Eq.(|2|), is of order unity per particle. The typical distance between ions scales as f ~ S^/^ (for 
S > 1) or like f ^ 2^'^ (for S < 1 and assuming a liquid-like structure) and thus the repulsive 
energy between two ions at this distance, as determined by the first term in Eq.(g), is of order 
^1/2 for S > 1 or S^/'^ for S < 1. It thus follows on rather general grounds that for a large 
coupling parameter, S > 1, the ionic structure is dominated by mutual repulsions, suggesting 
crystallization, while for weak coupling S < 1 these repulsions should be rather unimportant. 
In Fig.l we show ion-distribution snapshots from Monte-Carlo simulations for S = 0.1, 10, 10*. 
For small S, Fig. la, repulsive ion-ion interactions are indeed not playing a dominant role, the 
ion distribution is rather disordered and mean-field theory should work, since each ion moves 
in a diffuse cloud of neighboring ions. For large 2, on the other hand, ion-ion repulsions are 
strong and ion-ion distances are large compared to the distance from the wall. In Fig.lc the 
ions form a flat layer on the charged wall. A two-dimensional one-component plasma is known 
to crystallize for plasma parameter F « 125[^. From the definition of the two-dimensional 
plasma parameter llG] we obtain the relation S = 2r^. This leads to a crystallization threshold 
(in units of our coupling parameter) of S w 31000 meaning that all three systems in Fig.l are 
below the crystallization threshold. However, mean-field theory is expected to break down, at 
least for the system shown in Fig.lc, because each ion moves, though confined by its immediate 
neighbors in the lateral directions, almost independently from the other ions along the vertical 
direction (which constitutes the soft mode). Let us now see how these notions can be made 
concrete in a field-theoretic formulation. 

The partition function of a system of N counter-ions interacting via Coulomb interactions 
v(r) = 1/r among themselves and with a distribution cr(r) is given by 



N 
i 



— Jl / dr^ exp } -q^lB Y^ v{vj - Vk) + qts / drCT(r) ^ v{y ~ r^) + ^ h{Yj) 
■ i=i '' [ i<k •' j j 

'(3) 

where the field h has been added to calculate density distributions later on. We implicitly 
assume the configurational integral to span the upper half space {z > 0) only. At this point 
we employ a Hubbard-Stratonovitch transformation, similar to previous implementations of a 
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Fig. 1. - Snapshots of distributions of 75 ions for different values of the coupling constant H from 
Monte-Carlo simulations, a) Weak-coupling limit where Poisson-Boltzmann theory is accurate, H = 
0.1, b) intermediate coupling regime, H — 10, and c) strong-coupling limit, H — 10000. Note that in 
a) 8 particles which are located far away from the charged surface are not shown. 



field theory for charged system R, O], which leads to 

^"^ " M / ? ''''P {" 2^ / '^^^'■>W""'(^ - '■')<^('"') + *y drc7(r)0(r)| J 



-\N 



^j.Qh(r)-iq,j}{r) 

(4)" 

where we introduced the notation Zy = \/dct v. For simplicity, we perform a Legendre 
transformation to the grand-canonical ensemble, Q = '^^ ^^ ^n, and thereby introduce the 
fugacity A. The inverse Coulomb operator follows from Poisson's law as v~^{r) = — V^(5(r)/47r. 
The grand-canonical partition function can therefore be written as 



s = / -^ exp <; - 



dr 



[V0(r 



^l2 

^ - ia{v)(t,{v) - A6i(0)e''('')~'«'^('') 



(5) 



where 6{z) = 1 for z > and zero otherwise. Next we rescale the action, similarly to 
our rescaling analysis of the Hamiltonian in the beginning. All lengths are rescaled by the 
Gouy-Chapman length, r = /if, the fluctuating field is rescaled by the valency, (j) = cfs/q. We 
also use the explicit form a{v) = as6{z) for the fixed charge distribution. The result is 



S = / ^ exp 



1 



SttS 



df 



[V^(f)]^ - 4i(5(z)^(f) - 4A6i(z) 



5^e'i(r)— ' 



(6) 
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where the rescaled fugacity A is defined by 

A = 27rAM^S = -A^. (7) 

The expectation value of the counter-ion density, {p{r)), follows by taking a functional deriva- 
tive with respect to the generating field h, {p{r)) — S\nQ/dh{r)p^, giving rise to 



2TTiB<y^s 



= A(c-*'^^^-'). (8) 



The normalization condition for the counter- ion distribution, fi J dzp{z) = cFs/q, which follows 
directly from the definition of the grand-canonical partition function, leads to 

/•OO 

A / d5(e-*"^(^)) = 1. (9) 



"'0 

This is an important equation since it shows that the expectation values of the fugacity term in 
Eq. (pi) is bounded and of the order of unity per unit area. Let us first repeat the saddle-point 
analysis, which, because of the structure of the action in Eq.(^), should be valid for S <C 1. 
The saddle-point equation reads 

^0- = 2zAe-*^(^) (10) 

with the boundary condition d(f>{z)/dz ~ — 2i at z = 0. The solution of this differential 
equation is 



i^{S) = 2lii(l + A^^^s\ (11) 



while the boundary condition leads to A = 1 , which shows that the saddle-point approximation 
is indeed valid in the limit S <C 1. Combining Eqs.(|§) and (pi]), the density distribution of 
counter ions is given by the well-known result 



27r£Bf72 (1 + z)2 

Let us now consider the opposite limit, when the coupling constant S is large. In this case, 
the saddle-point approximation breaks down, since the prefactor in front of the action in 
Eq.(|6|) becomes small. This has been seen in a systematic loop-wise expansion around the 
saddle point, where the corrections to the saddle-point solution are proportional to powers 
of S and thus become largeH. However, from the field-theoretic partition function Eq.(0), 
it is self-evident what has to be done in this limit. Since the fugacity term is bounded, as 
evidenced by Eq.(P), one can expand the partition function (and also all expectation values) 
in powers of A/S. Upon Legendre transformation to the canonical ensemble, this gives the 
standard virial expansion. The normalization condition Eq. (m can be solved by an expansion 
of the fugacity as A = Aq + Ai/S + . . ., which leads to an expansion of the density profile 
with the small parameter 1/S. While the standard virial expansion fails for homogeneous bulk 
systems because of infra-red divergences, these divergences are renormalized for the present 
case of inhomogeneous distribution functions via the normalization condition Eq. (0) as we will 
now demonstrate. To leading order in this expansion, the rescaled density is 

(p(f)) _ , _ r s^,,^, , 1 



2ti(.bc^s 



= Aexp I -^v(O) + ^ / dn;(r - r)6{z) \ . (13) 
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Fig. 2. - Rescaled counter-ion density distribution p/2-K(.BUg as a function of the rescaled distance 
from the wall z/ p,. The inset shows Monte Carlo data for coupling constants H = 10^, 10**, 100, 10, 1, 0.1 
(from bottom to top) in a double- logarithmic plot, the main figure shows data for H = 0.1 (open 
diamonds), H = 10 (open triangles) and H = 10'* (open stars). The solid and broken lines denote the 
Poisson-Boltzmann and strong-coupling predictions, Eqs.(FL3) and (lla), respectively. All data were 
obtained with 75 particles and 10^ Monte Carlo steps (MC5), exceptthe data for H = 0.1 where 600 
particles were simulated. Error bars are smaller than the symbol size if not shown. 



From the normalization condition Eq.(p|) we obtain 

1 



Ao = exp I -v{Q) - ^ / drv{v)5{z) 

and thus the density distribution is to leading order given by 

(P(?)) _^-z 



(14) 



(15) 



The density profile in the strong coupling limit S — > oo is thus given by a simple exponential, 
while the profile in the weak-coupling (Poisson-Boltzmann) limit S -^ 0, Eq.dT2), is given by a 
power law. Both distributions give the same density at the wall, which is a trivial consequence 
of the contact- value theorem since the net force on the wall is zeroj£|. Just to avoid confusion 
at this stage, we stress that the exponential density profile Eq.([l5) has nothing to do with 
the Debye-Hiickel approximation. It is true that an exponential density profile (though with 
a different density contact value and thus violating the contact-value theorem) also follows 
from linearizing the differential equation Eq.(nG). However, the linearized solution can never 
be a more faithful representation of the true density profile than the full non-linear solution, 
and secondly, nothing in the saddle-point equation indicates that linearization should be valid 
for large values of S, because Eq.(nfl) does not depend on S. It is the saddle-point approach 
itself (the framework within which the Debye-Hiickel approximation can be formulated) which 
becomes invalid when S becomes large. 
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Fig. 3. - Results for the rescaled counter-ion density distribution pJ^-kIbo'I as a function of the 
rescaled distance from the wall zj [i for H — 0.1 and 5 particles and 10® MCS (open stars), 15 particles 
and lO'' MCS (open triangles), 35 particles and lO'' MCS (filled stars), and 600 particles and 10^ MCS 
(open diamonds) . Finite size are negligible except for the simulations using 5 particles and (to a lesser 
degree) 15 particles. 

In Fig. 2 we show counter-ion density profiles obtained using Monte Carlo simulations 
for various values of the coupling parameter S. As can be seen, the PB density profile 



Eq.(12) is only realized for S < 1, while the strong-coupling profile Eq.(|l5|) is indeed the 
asymptotic solution and agrees with simulation results for S > 10^. Experimentally, a coupling 
parameter S = 100, which is already quite close to the strong-coupling limit (see Fig. 2 inset), 
is reached with divalent ions for a surface charged density Og ~ 3.6nm~^, which is feasible 
with compressed charged monolayers, and with trivalent counter ions for <Ts w Inrn^'^, which 
is a typical value. The strong-coupling limit is therefore experimentally accessible. In our 
simulations we employ periodic boundary conditions in the lateral directions, while the ions are 
unconfined vertically. Forces and energies are calculated following the approach by Leknerjl2| 
and Sperbp3[ where the sum over the periodic images is transformed into a rapidly converging 
expansion in terms of Bessel functions ImI . Since the number of particles in the simulation, N, is 
related to the linear box size (in units of the Gouy-Chapman length fi) by L = L/ fi = \/2ttNE., 
effects due to the finite number of particles and effects due to the finite lateral box size are 
connected. In Fig. 3 we show MC density profiles for a fixed coupling constant S = 0.1 and for 
various number of particles. For all but the smallest systems with 5 and 15 particles finite-size 
effects are negligible. 

In summary, we have derived within a field-theoretic framework the strong-coupling theory 
for counter-ion distributions, valid in the limit of large valencies and/or low temperatures. It 
corresponds to a standard virial expansion of the field-theoretic action, which gives meaningful 
results because long-range divergences, which spoil the free-energy virial expansion for bulk 
systems and lead to non-analytic terms [[15|, are subtracted by a fugacity renormalization. This 
theory thereby complements the Poisson-Boltzmann approach, which works for the oppposite 
limit of weak coupling. We note that an exponential density profile (albeit with a different 
prefactor) has been derived in the strong-coupling limit by Shklovskii using a heuristic model 
where ions bound to the wall are in chemical equilibrium with free ions[pl. 

Next-leading order results for the density profile within our strong-coupling theory, and 
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comparison of these corrections with Monte-Carlo simulations, will be presented in a future 
publication. Also, generalization to other geometries, including the pressure between two 
charged walls, the presence of salt, dielectric boundaries, and hard-core interactions are 
straightforward within our theory and will be published soon. For the case of a cylinder 
with a radius R and linear charge density t we obtain that the strong-coupling theory is valid 
for Sc 3> 1 with the cylindrical coupling parameter defined as Sc = q^(.%T j R. The Manning 
condensation threshold is reflected by the fact that the density distribution is not normalizable 
for tIb < 1- For t£b > 1 only part of the counter ions are condensed and the normalization 
condition Eq.(^) is modified. The result for the density profile in the limit S^ ^ 1 is 

and thus exhibits purely algebraic decay, in contrast to the PB prediction [p6[ but in agreement 
with Shklovskii's results[p|. Finally, we note that the strong-coupling limit not necessarily 
entails unrealistically small Gouy-Chapman lengths. Writing the Gouy Chapman length as 
M = ^-B9^/2, it is clear that /j, can be fixed at a moderate value for highly valent ions even when 
S 3> 1. In this case, however, one would have to worry about the effect of discrete surface 
charges. 
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